sa.nps=function(y, x, xnet, tau=0.5, start, sigma, delta, rt=0.15, eps=1.0E-6,
                   nt=20, ns=10, int=T, neps=4, maxevl=10000000, maxhill=1000,
                   sl, vm, seed1=10001, seed2=10002, t0, lbc, ubc,
                   sgn=1, sa=T, nr=T, hill=T, bc=T, rho=1)
{
  
  #**********************************************************************
  #				PART I
  #   Check the inputs of the function, read data and set parameters.
  #**********************************************************************
  if(!is.loaded("npsms")) dyn.load("npsms.dll");
  is.loaded("npsms")
  
  tau = as.vector(tau)					# Read the quantile(s) to be estimated.
  p = length(tau)					# Number of quantiles  to be estimated.
  y = as.vector(y)           				# Read the dependent latent (0,1) variable.
  x = as.matrix(cbind(xnet,x))           				# Read the expalnatory variables.
  if(int) x = cbind(1,x)     				# Add an intercept to the model.
  nobs = nrow(x)             				# Number of observations.
  k = ncol(x)                				# Number of explanatory variables.
  if(int) {if(k<3) stop()}				# Check that there is a regressor to standardize.
  else {if(k<2) stop()} 
  
  
  # Get the names of the variables.
  if(int) xnames <- c("intercept","network", dimnames(x) [[2]] [3:k])
  else xnames <- c("network",dimnames(x) [[2]] [2:k])
  
  
  if(missing(sigma)) {sigma = rep(nobs^(-1/5),p)} 		# Set the smoothing parameter.
  else if(length(sigma)!=p) stop("sigma and tau should have the same lenght.")
  if(missing(delta)) {delta = rep(0.1,p)}                	# Set delta.
  else if(length(delta)!=p) stop("delta and tau should have the same lenght.")
  sigmad = sigma^(delta)         			# Set the smoothing parameter for Bias Estimation.
  
  lb = rep(lbc,k-1)					# Set lower bound of the parameter space.
  lb[2]=0;
  ub = rep(ubc,k-1)					# Set upper bound of the parameter space.
  
  if(missing(start)) {start = matrix(rep(rep(1,k),p),k,p)} 			# Set starting values.
  else {start = as.matrix(start)   
  if(ncol(start)!=p) stop("number of columns of start should equal the lenght of tau.")
  if(nrow(start)!=k) stop("number of rows of start should equal the number of columns of x.")
  }
  #**************************************************
  # Check that the starting values are within bounds.
  #**************************************************
  for(i in 1:(k-1)) {
    for(m in 1:p){
      if(start[i,m]<lb[i] | start[i,m]>ub[i]) stop(
        cat("Initial values (start) out of bounds (lb, ub).","\n")) }}
  
  if(missing(sl)) {sl = rep(2, k-1)}
  if(missing(vm)) {vm = rep(1, k-1)}
  
  #************************************************
  # Check the seed of the random number generator.*
  #************************************************
  if(seed1<0 | seed1>31328 | seed2<0 | seed2>30081)
    stop("The first number seed (seed1) must have a value between 0
         and 31328 and the second value seed (seed2) must have a value between
         0 and 30081.")
  
  #*****************************
  # Set the initial temperature.
  #*****************************
  if(missing(t0)) {
    f = rep(sum((y-(1-0.5))*
                  pnorm((as.matrix(x[,1:(k-1)],nobs,k-1)%*%(runif(k-1,lb[[1]],ub[[1]]))+sgn*x[,k])/sigma[1])),
            50000)
    t0 = as.double(sqrt(var(f)))
  }
  else {if(t0<0) stop(" The initial temperature (t0) is not positive. 
                      Reset t0 to a positive number, or let the routine pick the initial 
                      temperature automatically.")}
  
  
  

  
  #**********************************************************************
  #				PART II                               *
  #			Maximize the objective function.              *
  #**********************************************************************
  
  betavec = NULL
  objvec = NULL
  
  for(m in 1:p) {
    taum = tau[m]
    #print(m)
    #***********************************************
    # Evaluate the objective at the starting values.
    #***********************************************
    
    fstart = 1/(nobs)*sum((y-(1-taum))*
                            pnorm((as.matrix(x[,1:(k-1)],nobs,k-1)%*%start[1:(k-1),m]+sgn*x[,k])/sigma[m]))
    # cat("\n",
    #     "*******************************************************************"     , "\n",
    #     "            Estimation of the", taum, "Quantile Model"                   , "\n",
    #     "*******************************************************************"     , "\n",
    #     "                 ", xnames, "Objective"                                  , "\n",
    #     "-------------------------------------------------------------------"     , "\n",
    #     " Initial Values  ", format(round(start[,m], 6)), format(round(fstart,6)) , "\n"
    # )
    
    beta = start[,m]
    fmax = fstart
    
    if(sa){
      #**********************************************************************
      # 		Invoke the SA algorithm for function optimization.    *
      #**********************************************************************
      # cat("\n",
      #     " Invoking the SA maximization routine."                                 , "\n",
      #     " Please wait..."                                                        , "\n"
      # )
      
      out2 = .Fortran("sa",
                      as.double(y),
                      as.double(x),
                      as.integer(nobs),
                      as.integer(k),
                      as.double(start[,m]),
                      taum = as.double(taum),
                      as.logical(T),
                      as.double(rt),
                      as.double(eps),
                      as.integer(ns),
                      as.integer(nt),
                      as.integer(neps),
                      as.integer(maxevl),
                      as.double(lb),
                      as.double(ub),
                      as.double(sl),
                      as.double(seed1),
                      as.double(seed2),
                      as.double(t0),
                      as.double(vm),
                      coef = as.double(start[,m]),
                      obj = as.double(0),
                      as.integer(0),
                      as.integer(0),
                      as.integer(0),
                      as.integer(0),
                      as.double(start[,m]),
                      as.double(start[,m]),
                      as.integer(start[,m]),
                      aflag = as.integer(0),
                      sigmam = as.double(sigma[m]),
                      as.double(sgn)
      )
      
      if(out2$aflag!=0)
        warning("Too many function evaluations. Consider increasing maxevl
                or eps, or decreasing nt or rt. These results are likely
                to be poor.")
      
      # cat("\n", " SA finished", "\n")
      
      beta = as.matrix(out2$coef)
      fmax = 1/(nobs)*sum((y-(1-taum))*
                            pnorm((as.matrix(x[,1:(k-1)],nobs,k-1)%*%beta[1:(k-1)]+sgn*x[,k])/sigma[m]))
      
      # cat("\n",
      #     "Estimates After SA"                              , "\n",
      #     "------------------"                              , "\n",
      #     xnames, "Objective"                               , "\n",
      #     format(round(beta, 8)), format(round(fmax,8))     , "\n"
      # )
    }								# end if for SA
  	    
    
    
    betavec = cbind(betavec,beta)
    objvec = cbind(objvec,fmax)
  }							# End of tau iteration.
 
  
  list(coef=betavec,sigma=sigma)
}